/*
 *  VortMesh.cpp
 *  VortPan
 *
 *  Created by Gregor Veble on 6/29/10.
 *  Copyright 2010  . All rights reserved.
 *
 */

#include <assert.h>
#include <algorithm>
#include <fstream>
#include <sstream>
#include <iostream>
#include "VortMesh.h"

using namespace std;
using namespace VortPan;

VortMesh::VortMesh() : Baked(false), NSolids(0) {
}

VortMesh::VortMesh(const VortMesh& Other){
	*this=Other;
}


VortMesh::~VortMesh(){
}

VortPointReference VortMesh::AddPoint(const VortPoint& P){
	PointStorage.push_back(P);

	return GetPointReference(PointStorage.size()-1);
}

int VortMesh::AddPointPolygon(VortPolygon& P){
	int res=0;
	for (PointReferenceVector::iterator it=P.PointList.begin(); it!=P.PointList.end(); it++) {
		if (it->Index<0) { // Warning; this is for polygons that are only defined by pointers to points; better not use
			int FoundPointIndex=GetPointIndex((*it));
			if(FoundPointIndex>=0){
				it->Index=FoundPointIndex;
			} else {
				res=-1;
			}
		}
	}
	if(res==0){
		PolygonStorage.push_back(P);
		return PolygonStorage.size()-1;
	} else {
		return res;
	}
}

VortPointReference VortMesh::GetPointReference(int index){
	VortPointReference rPoint(PointStorage[index],index);
	return rPoint;
}

int VortMesh::Bake(){
	RenewAllPointReferences();
	BuildEdges();
//	cout<<"PointSize="<<PointStorage.size()<<endl;
	Baked=true;
	return 0;
}

int VortMesh::GetPointIndex(const VortPointReference& P){

	int Size=PointStorage.size();
	
	for (int i=0; i<Size; i++) {
		if (P.pPoint==&(PointStorage[i])) {
			return i;
		}
	}
	
	return -1;
}

void VortMesh::Clear(){
	Baked=false;
	PolygonStorage.clear();
	PointStorage.clear();
}

Vector3D RandDir={0.810275, 0.159144, -0.564029};
	
struct WeldHelper {
	int OldPointIndex;
	int EqualTo;
	int NewPointIndex;
	re Coord;
};

bool WeldComparison(const WeldHelper& W1, const WeldHelper& W2){
	return (W1.Coord<W2.Coord);
}

bool OriginalOrderComparison(const WeldHelper& W1, const WeldHelper& W2){
	return (W1.OldPointIndex<W2.OldPointIndex);
}


int VortMesh::Weld(re Eps){
	
	
	PointVector WeldedPointStorage;
	
	vector<WeldHelper>  WeldArray;
	vector<WeldHelper>::iterator it;
	re Eps2=Eps*Eps;
	WeldHelper WH;
	
	WH.OldPointIndex=-1;
	WH.EqualTo=-1;
	WH.NewPointIndex=-1;
	WH.Coord=0;
	
	VortMesh WeldedMesh;
	
	WeldArray.assign(PointStorage.size(),WH);
	
	
	int i=0;
	for (it=WeldArray.begin(); it!=WeldArray.end(); it++) {
		it->OldPointIndex=i;
		it->EqualTo=i;
		it->NewPointIndex=10000;
		it->Coord=PointStorage[i].Position*RandDir;
		i++;
	}
	
	sort(WeldArray.begin(), WeldArray.end(), WeldComparison);
	
	for (it=WeldArray.begin(); it!=WeldArray.end(); it++) {
		VortPoint& P1=PointStorage[it->OldPointIndex];
		for(vector<WeldHelper>::iterator it2=it+1;(it2!=WeldArray.end()) && ((it2->Coord-it->Coord)<=Eps) ;it2++){
			VortPoint& P2=PointStorage[it2->OldPointIndex];
			if (Square(P2.Position-P1.Position)<=Eps2) {
				it2->EqualTo=it->EqualTo;
			}
		}
	}
	
	sort(WeldArray.begin(), WeldArray.end(), OriginalOrderComparison);

	i=0;
	for (it=WeldArray.begin(); it!=WeldArray.end(); it++) {
		if (it->EqualTo==it->OldPointIndex) {
			it->NewPointIndex=i;
			WeldedMesh.AddPoint(PointStorage[it->OldPointIndex]);
			i++;
		}
	}
	
	
	int NWelded=0;
	for (it=WeldArray.begin(); it!=WeldArray.end(); it++) {
		if (it->EqualTo!=it->OldPointIndex) {
			it->NewPointIndex=WeldArray[it->EqualTo].NewPointIndex;
			NWelded++;
		} 
	}
	
	
	for (PolygonVector::iterator it1=PolygonStorage.begin(); it1!=PolygonStorage.end(); it1++) {
		VortPolygon& rOldP=*it1;
		VortPolygon P(rOldP.WakePoly);
		
		P.SetCenterPoint(rOldP.CenterPoint);

		for (PointReferenceVector::iterator itref=rOldP.PointList.begin(); itref!=rOldP.PointList.end(); itref++) {
			int NewInd=WeldArray[itref->Index].NewPointIndex;
			P.AddPoint(WeldedMesh.GetPointReference(NewInd));
		}
		WeldedMesh.AddPointPolygon(P);
	}
	
	*this=WeldedMesh;
	
	Bake();
	
//	cout<<"NWelded = "<<NWelded<<endl;
	return NWelded;
}

void VortMesh::AddSimplePolygon(vector<Vector3D>& CornerList){
	bool OldBaked=Baked;
	
	Baked=false;
	int index=PointStorage.size();
	
	VortPolygon Poly;
	VortPoint Point;
	for (vector<Vector3D>::iterator it=CornerList.begin(); it!=CornerList.end(); it++) {
		Point.Position=*it;
		AddPoint(Point);
		Poly.AddPoint(GetPointReference(index));
		index++;
	}
	
	Baked=OldBaked;
}

void VortMesh::FinishSimple(){
	Bake();
}

void VortMesh::RenewAllPointReferences(){
	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		PointReferenceVector& PL=it->PointList;
		for (PointReferenceVector::iterator itr=PL.begin(); itr!=PL.end(); itr++) {
			itr->pPoint=&(PointStorage[itr->Index]);
		}
	}
}

void VortMesh::BuildEdges(){
	EdgeStorage.clear();
	
	int PolyIndex=0;
	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		for (PointReferenceVector::iterator itp=it->PointList.begin(); itp!=it->PointList.end(); itp++) {
			PointReferenceVector::iterator itp1=itp+1;
			if (itp1==it->PointList.end()) {
				itp1=it->PointList.begin();
			}
			VortEdge E(*itp,*itp1);
			
			PolygonReferenceVector& PRV=EdgeStorage[E].PolygonReferences; 
			//this also creates EdgeStorage[E] if it doesn't exist yet
			
			EdgeStorage[E].Problematic=false;
			
			VortPolygonReference PolyRef(*it,PolyIndex);
			PolyRef.EdgeDirection=E.Direction;
			
			PRV.push_back(PolyRef);
			
		}
		PolyIndex++;
	}
	
//	cout<<"NEdges = "<<EdgeStorage.size()<<endl;
}

int VortMesh::CheckEdges(bool Symmetric, re Eps, SymmetryPlaneVector SymmetryPlanes){
	int NWrong=0;
	
	int iptemp=0;
	
	for (EdgeMap::iterator it=EdgeStorage.begin(); it!=EdgeStorage.end(); it++){
		PolygonReferenceVector& PRV=it->second.PolygonReferences;
		it->second.Problematic=false;
		int NPoly=PRV.size();
		int NWake=0;
		int iWake=0;
		iptemp++;
		switch (NPoly) {
			case 1:
				if (!PRV[0].pPolygon->WakePoly) {
					if (Symmetric) {
						bool TestedPlanes=false;
						for (SymmetryPlaneVector::iterator its=SymmetryPlanes.begin(); its!=SymmetryPlanes.end(); its++) {
							VortSymmetryPlane& Plane=*its;
							re Off1=Plane.PlaneDistance(it->first.PointRef1.pPoint->Position);
							re Off2=Plane.PlaneDistance(it->first.PointRef2.pPoint->Position);
							if ((rabs(Off1)<Eps)&&(rabs(Off2)<Eps)) {
								TestedPlanes=true;
							}
						}
						if (!TestedPlanes) {
							NWrong++;
							it->second.Problematic=true;
						}
					} else {
						NWrong++;
						it->second.Problematic=true;
					}

				}
				break;
			case 2:
				if (PRV[0].EdgeDirection==PRV[1].EdgeDirection) {
					NWrong++;
					it->second.Problematic=true;
				} else{
					if (PRV[0].pPolygon->WakePoly!=PRV[1].pPolygon->WakePoly) {
						NWrong++;
						it->second.Problematic=true;
					}
				}
				break;
			case 3:
				for (int i=0; i<3; i++) {
					if (PRV[i].pPolygon->WakePoly) {
						NWake++;
						iWake=i;
					}
				}
				if (NWake!=1) {
					NWrong++;
					it->second.Problematic=true;
				}
				if ((PRV[(iWake+1)%3].pPolygon->WakePoly)||(PRV[(iWake+2)%3].pPolygon->WakePoly)) {
					NWrong++;
					it->second.Problematic=true;
				} else {
					if (PRV[(iWake+1)%3].EdgeDirection==PRV[(iWake+2)%3].EdgeDirection) {
						NWrong++;
						it->second.Problematic=true;
					}
				}
				break;

			default:
				NWrong++;
				it->second.Problematic=true;
				break;
		}
	}

	return -NWrong;
}


int VortMesh::LabelSolids(){

	NSolids=0;
	
	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		it->pSibling=&*it;
		it->SolidIndex=-1;
	}

	for (EdgeMap::iterator it=EdgeStorage.begin(); it!=EdgeStorage.end(); it++){
		PolygonReferenceVector& PRV=it->second.PolygonReferences;
		int NPoly=PRV.size();
		for (int i=1; i<NPoly; i++) {
			VortPolygon* pTop=PRV[i].pPolygon->GetTopmostSibling();
			pTop->pSibling=PRV[0].pPolygon->GetTopmostSibling();
		}
		
	}
	
	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		VortPolygon* pTop=it->GetTopmostSibling();
		if(pTop->SolidIndex==-1){
			pTop->SolidIndex=NSolids;
			NSolids++;
		}
		it->SolidIndex=pTop->SolidIndex;
	}
	
	
	return NSolids;
}

int VortMesh::AddMesh(VortMesh& Other){
	if(&Other==this){
		return -1;
	}
	
	RenewAllPointReferences();
	Other.RenewAllPointReferences();
	
	int Offset=PointStorage.size();

	PointStorage.insert(PointStorage.end(),Other.PointStorage.begin(),Other.PointStorage.end());
	
	for (PolygonVector::iterator it=Other.PolygonStorage.begin(); it!=Other.PolygonStorage.end(); it++) {
		for (PointReferenceVector::iterator itp=it->PointList.begin(); itp!=it->PointList.end(); itp++) {
			itp->Index+=Offset;
		}
	}
	
	PolygonStorage.insert(PolygonStorage.end(),Other.PolygonStorage.begin(),Other.PolygonStorage.end());

	Bake();
	
	return 0;
}

int VortMesh::ImportRawData(std::string FName, bool WakeMesh,  re ScaleFactor, int IndOffset, bool OnlyPos){
	Clear();
	
	ifstream inputfile(FName.c_str());
	
	if(!inputfile.is_open()){
		return -1;
	}
	
	int NPoints;
	
	inputfile>>NPoints;

	VortPoint P;
	
	for(int i=0; i<NPoints; i++){
		inputfile>>P.Position.v[0];
		inputfile>>P.Position.v[1];
		inputfile>>P.Position.v[2];
		if(!OnlyPos){
            inputfile>>P.Normal.v[0];
            inputfile>>P.Normal.v[1];
            inputfile>>P.Normal.v[2];
            inputfile>>P.UV.v[0];
            inputfile>>P.UV.v[1];
        }		
		P.Position=ScaleFactor*P.Position;
		
		AddPoint(P);

	}

	int NPolys;
	
	inputfile>>NPolys;
	
	for(int i=0; i<NPolys; i++){
		int NVert;
		
		inputfile>>NVert;
		
		VortPolygon P(WakeMesh);
		VortPoint Dummy;
		for (int k=0; k<NVert; k++) {
			int PointInd;
			inputfile>>PointInd;
            PointInd-=IndOffset;
			VortPointReference VP(Dummy,PointInd);
			P.AddPoint(VP);
		}
		P.EliminateDegenerateEdges();
		AddPointPolygon(P);
		
	}
	
	Bake();

	RecalculateCenters();
	
	return 0;
}



int VortMesh::LoadPLY(std::string FName, bool WakeMesh, re ScaleFactor){
	
	Clear();
	
	int xoff=0,yoff=0,zoff=0,faceoff=0,vertind=0,faceind=0;
	int elementindex=0;
	string str;
	int numvert,numface;
	ifstream inputfile(FName.c_str());

	if(!inputfile.is_open()){
		return -1;
	}
	
	string key;
	bool ReadingVert=false;
	bool ReadingFace=false;
	do {
		getline(inputfile,str);
		stringstream sstr(str);
		
		
		sstr>>key;
		if(key=="element"){
			sstr>>key;
			if((key=="vertex")){
				if(elementindex!=0){
					return -6;
				}
				ReadingVert=true;
				ReadingFace=false;
				sstr>>numvert;
			} else {
				if(key=="face"){
					if(elementindex!=1){
						return -7;
					}
					ReadingVert=false;
					ReadingFace=true;
					sstr>>numface;
				} else {
					if(elementindex<2){
						return -3;
					}
				}
			}
			elementindex++;
		}
		
		if(key=="property"){
			sstr>>key;
			if(ReadingVert){
				sstr>>key;
				if(key=="x"){
					xoff=vertind;
				}
				if(key=="y"){
					yoff=vertind;
				}
				if(key=="z"){
					zoff=vertind;
				}
				vertind++;
			}
			if(ReadingFace){
				if(key=="list"){
					sstr>>key;
					sstr>>key;
					sstr>>key;
					if(key=="vertex_index"){
						faceoff=faceind;
					}
				}
				faceind++;
			}
		}

	} while ((!inputfile.eof())&&(key!="end_header"));

	if(inputfile.eof()){
		return -2;
	}
	
	int i;
	for (i=0; (i<numvert)&&(!inputfile.eof()); i++) {
		getline(inputfile,str);	
		stringstream sstr(str);
		VortPoint P;
		for (int j=0; j<vertind; j++) {
			if(j==xoff){
				sstr>>P.Position.v[0];
			} else {
				if(j==yoff){
					sstr>>P.Position.v[1];
				} else {
					if(j==zoff){
						sstr>>P.Position.v[2];
					} else {
						sstr>>str;
					}
				}
			}
		}
		P.Position=ScaleFactor*P.Position;
		AddPoint(P);
	}

	if(i!=numvert){
		Clear();
		return -4;
	}
	
	
	for (i=0; (i<numface)&&(!inputfile.eof()); i++) {
		getline(inputfile,str);			
		stringstream sstr(str);
		for (int j=0; j<faceind; j++) {
			if (j==faceoff) { 
				int nv;
				sstr>>nv;
				VortPolygon P(WakeMesh);
				VortPoint Dummy;
				for (int k=0; k<nv; k++) {
					int PointInd;
					sstr>>PointInd;
					VortPointReference VP(Dummy,PointInd);
					P.AddPoint(VP);
				}
				P.EliminateDegenerateEdges();
				AddPointPolygon(P);
			} else {
				sstr>>str;
			}
		}		
	}
	if (i!=numface) {
		return -5;
	}
	
	Bake();
	
	RecalculateCenters();
	

	return 0;
}

int VortMesh::SavePLY(std::string FName){
	ofstream outputfile(FName.c_str());
	
	if(!outputfile.is_open()){
		return -1;
	}
	
	outputfile<<"ply"<<endl;
	outputfile<<"format ascii 1.0"<<endl;
	outputfile<<"comment created by VortPan, Pipistrel d.o.o. Ajdovscina, www.pipistrel.si"<<endl;
	outputfile<<"element vertex "<<PointStorage.size()<<endl;
	outputfile<<"property float32 x"<<endl;
	outputfile<<"property float32 y"<<endl;
	outputfile<<"property float32 z"<<endl;
	outputfile<<"element face "<<PolygonStorage.size()<<endl;
	outputfile<<"property list uint8 uint32 vertex_indices"<<endl;
	outputfile<<"end_header"<<endl;

	Vector3D Point;
	for (PointVector::iterator it=PointStorage.begin(); it!=PointStorage.end(); it++) {
		Point=it->Position;
		outputfile<<Point.v[0]<<" "<<Point.v[1]<<" "<<Point.v[2]<<endl;
	}
	
	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		outputfile<<it->PointList.size();
		for (PointReferenceVector::iterator itp=it->PointList.begin(); itp!=it->PointList.end(); itp++) {
			outputfile<<" "<<itp->Index;
		}
		outputfile<<endl;
	}
	
	outputfile<<endl;
	
	return 0;
}


void VortMesh::RecalculateCenters(re Eps){
	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		it->RecalculateAverageCenterPoint(Eps);
	}
}

int VortMesh::ExportCenterPoints(std::string FName){
	ofstream of(FName.c_str());
	
	if(!of.is_open()){
		return -1;
	}
	
	fixed(of);

	for (PolygonVector::iterator it=PolygonStorage.begin(); it!=PolygonStorage.end(); it++) {
		of<<it->GetCenterPosition().v[0]<<" "<<it->GetCenterPosition().v[1]<<" "<<it->GetCenterPosition().v[2]<<endl;
	}
	
	return 0;
}


int VortMesh::ExportMathematicaList(std::string FName){
	ofstream of(FName.c_str());
	
	if(!of.is_open()){
		return -1;
	}
	
	fixed(of);
	
	
	of<<"{"<<endl;
	for (int k=0; k<PolygonStorage.size(); k++) {
		of<<"{";
		for (int i=0; i<PolygonStorage[k].GetNPoints(); i++) {
			VortPoint& VP=PolygonStorage[k].GetPoint(i);
			of<<"{";
			for (int j=0; j<3; j++) {
				of<<VP.Position.v[j];
				if (j!=2) {
					of<<", ";
				}
			}
			of<<"}";
			if (i<(PolygonStorage[k].GetNPoints()-1)) {
				of<<", ";
			}
		}
		of<<"}";
		if(k<(PolygonStorage.size()-1)){
			of<<", ";
		}
		of<<endl;
	}
	of<<"}"<<endl;
	
	return 0;
}

int VortMesh::ExportBadEdges(std::string FName){
	ofstream of(FName.c_str());
	
	if(!of.is_open()){
		return -1;
	}
	
	fixed(of);

	int NBad=0;
	for (EdgeMap::iterator it=EdgeStorage.begin(); it!=EdgeStorage.end(); it++) {
		if(it->second.Problematic){
			NBad++;
			Vector3D Pos1=it->first.PointRef1.pPoint->Position;
			Vector3D Pos2=it->first.PointRef2.pPoint->Position;
			of<<Pos1.v[0]<<" "<<Pos1.v[1]<<" "<<Pos1.v[2]<<" ";
			of<<Pos2.v[0]<<" "<<Pos2.v[1]<<" "<<Pos2.v[2]<<endl;
		}
	}
	
	return NBad;
}

int VortMesh::Mirror(VortSymmetryPlane Plane){
	
	for (PointVector::iterator it=PointStorage.begin(); it!=PointStorage.end(); it++) {
		VortPoint& P=*it;
		P.Position=Plane.Reflect(P.Position);
	}
	
	
	for (PolygonVector::iterator it1=PolygonStorage.begin(); it1!=PolygonStorage.end(); it1++) {
		VortPolygon& P=*it1;

		P.CenterPoint.Position=Plane.Reflect(P.CenterPoint.Position);
		
		P.Normal=Plane. Mirror(P.Normal);

		reverse(P.PointList.begin(), P.PointList.end());
	}
	
	void BuildEdges();
	
	
	return 0;
	
}

VortMesh& VortMesh::operator=(const VortMesh& Mesh){
	if (&Mesh!=this) {
		PolygonStorage=Mesh.PolygonStorage;
		PointStorage=Mesh.PointStorage;
		EdgeStorage=Mesh.EdgeStorage;
		
		NSolids=Mesh.NSolids;

		RenewAllPointReferences();
		if(Mesh.Baked){
			Bake();
		}		
	}
	
	return *this;	
}


void VortMesh::CalculateMeshDerivative(std::vector<re>& Field, std::vector<Vector3D>& Derivative, re Parallel){
	Derivative.assign(PolygonStorage.size(),V0());
	assert(Field.size()>=PolygonStorage.size());
	re Val;
	
	for (int i=0; i<PolygonStorage.size(); i++) {
		VortPolygon& P=PolygonStorage[i];
		
		Derivative[i]=V0();
		
		if(!P.WakePoly){
			Vector3D D1=e1(),D2=e2();
			Vector3D& N=P.GetNormal();
			
			TriVectors(&N, &D1, &D2);
			
			
			
			for (PointReferenceVector::iterator itp=P.PointList.begin(); itp!=P.PointList.end(); itp++) {
				PointReferenceVector::iterator itp1=itp+1;
				if (itp1==P.PointList.end()) {
					itp1=P.PointList.begin();
				}
				VortEdge E(*itp,*itp1);	
				
				assert(EdgeStorage.find(E)!=EdgeStorage.end());
				
				VortEdgeInfo& EI=EdgeStorage[E];
				

				Val=0;//Field[i];
		
				if(EI.PolygonReferences.size()==2){
					int OtherI=-1;
					
					bool found=false;
					for (int j=0; j<2; j++) {
						if (EI.PolygonReferences[j].Index!=i) {
							OtherI=EI.PolygonReferences[j].Index;
						} else {
							found=true;
						}
					}
					assert(found);
					
					if (N*PolygonStorage[OtherI].GetNormal()>Parallel) {
						Val=HALF*(Field[OtherI]-Field[i]);
					}
				}
				
				Derivative[i]=Derivative[i]+Val*Vectorial(itp1->pPoint->Position-itp->pPoint->Position,N);
				
				
			}
						
			
			Derivative[i]=Derivative[i]/P.GetArea();
		}
	}
}


void VortMesh::CalculateMeshDerivativeLMS(std::vector<re>& Field, std::vector<Vector3D>& Derivative){
	Derivative.assign(PolygonStorage.size(),V0());
	assert(Field.size()>=PolygonStorage.size());

	for (int i=0; i<PolygonStorage.size(); i++) {
		VortPolygon& P=PolygonStorage[i];
		
		Derivative[i]=V0();
		
		if(!P.WakePoly){
			Vector3D D1=e1(),D2=e2();
			
			TriVectors(&P.GetNormal(), &D1, &D2);
			
			re a=0,b=0,d=0,e=0,f=0,x=0,y=0;
			
			
			for (PointReferenceVector::iterator itp=P.PointList.begin(); itp!=P.PointList.end(); itp++) {
				PointReferenceVector::iterator itp1=itp+1;
				if (itp1==P.PointList.end()) {
					itp1=P.PointList.begin();
				}
				VortEdge E(*itp,*itp1);	
				
				assert(EdgeStorage.find(E)!=EdgeStorage.end());
				
				VortEdgeInfo& EI=EdgeStorage[E];
				
				re ValDiff=0;
				
				Vector3D MidPointDiff;
				
				re Area=0;
				
				bool good=true;
				
				if((EI.PolygonReferences.size()==2)&&(EI.PolygonReferences[0].pPolygon->Normal*EI.PolygonReferences[1].pPolygon->Normal>0.5)){
					int OtherI=-1;
					
					bool found=false;
					for (int j=0; j<2; j++) {
						if (EI.PolygonReferences[j].Index!=i) {
							OtherI=EI.PolygonReferences[j].Index;
						} else {
							found=true;
						}
					}
					assert(found);
					
					ValDiff=(Field[OtherI]-Field[i]);
					MidPointDiff=PolygonStorage[OtherI].GetCenterPosition()-P.CenterPoint.Position;
					//MidPointDiff=TWO*(HALF*(itp1->pPoint->Position+itp->pPoint->Position)-P.CenterPoint.Position);
					Area=PolygonStorage[OtherI].GetArea();
				} else {
					MidPointDiff=TWO*(HALF*(itp1->pPoint->Position+itp->pPoint->Position)-P.CenterPoint.Position);
					ValDiff=0;
					Area=PolygonStorage[i].GetArea();
					good=false;
				}
				
				re MD2=MidPointDiff*MidPointDiff;
				
				re r1=MidPointDiff*D1;
				re r2=MidPointDiff*D2;
                
 				
				if(good){
					a+=r1*r1*Area/MD2;
					b+=r1*r2*Area/MD2;
					d+=r2*r2*Area/MD2;
					e+=r1*ValDiff*Area/MD2;
					f+=r2*ValDiff*Area/MD2;
				}
				
				
			}
			
			re Det=a*d-b*b;
			
			if(Det!=0){
				Det=ONE/Det;
				x=Det*(d*e-b*f);
				y=Det*(a*f-b*e);
			} 
			
			
			Derivative[i]=x*D1+y*D2;
		}
	}
}

void VortMesh::ReplaceAllWakePolys(Vector3D Direction, re StartSegment, re EndDistance, re Factor, re Eps){
	for (EdgeMap::iterator it=EdgeStorage.begin(); it!=EdgeStorage.end(); it++) {
		ReplaceWakePoly(it->first,it->second, Direction, StartSegment,EndDistance,Factor);
	}
	
	Bake();
	
	Weld(Eps);
}




void VortMesh::ReplaceWakePoly(VortEdge Edge, VortEdgeInfo& EdgeInfo, Vector3D Direction, 
							   re StartSegment, re EndDistance, re Factor){
	

	if (EdgeInfo.PolygonReferences.size()==3) {
		for (int i=0; i<EdgeInfo.PolygonReferences.size(); i++) {
			int index=-1;
			if (EdgeInfo.PolygonReferences[i].pPolygon->WakePoly) {
				index=i;
			}
			if (index>=0) {
				Direction=Normalize(Direction);
				VortPolygon& P=*EdgeInfo.PolygonReferences[index].pPolygon;
				int NPoint=4;
				if (P.PointList.size()==NPoint) {
					VortPointReference StartRef=Edge.PointRef1;
					if (EdgeInfo.PolygonReferences[index].EdgeDirection) {
						StartRef=Edge.PointRef2;
					} 
					int StartIndex=P.GetIndex(StartRef);
					int BeforeIndex=(StartIndex+NPoint-1)%NPoint;
					//StartIndex=(StartIndex+1)%P.PointList.size();
					
					re Distance=StartSegment;
					
					bool First=true;
					
					while (Distance<EndDistance) {

						VortPoint P1,P2;
						P1.Position=P.PointList[StartIndex].pPoint->Position+StartSegment*Direction;
						P2.Position=P.PointList[BeforeIndex].pPoint->Position+StartSegment*Direction;
						
						if (!First) {
							VortPointReference VPR1=AddPoint(P1);
							VortPointReference VPR2=AddPoint(P2); 
						
							P.AddPoint(VPR2,BeforeIndex);
							P.AddPoint(VPR1,BeforeIndex);
							StartIndex++;
							BeforeIndex++;
						} else {
							PointStorage[P.PointList[(StartIndex+1)%P.PointList.size()].Index].Position=P1.Position;
							PointStorage[P.PointList[(StartIndex+2)%P.PointList.size()].Index].Position=P2.Position;
							
							First=false;
							BeforeIndex+=P.PointList.size()-1;
							BeforeIndex%=P.PointList.size();
							StartIndex++;
											
						}
						StartSegment*=Factor;
						Distance+=StartSegment;
					}
				}
			}
		}
	}
}
